In this paper, the laminar incompressible flow equations are solved by an upwind LEAST-SQUARES MESHLESS METHOD. Due to the difficulties in generating quality meshes, particularly in complex geometries, a MESHLESS METHOD is increasingly used as a new numerical tool. The MESHLESS METHODs only use clouds of nodes to influence the domain of every node. Thus, they do not require the nodes to be connected to form a mesh and decrease the difficulty of meshing, particularly around complex geometries. In the literature, it has been shown that the generation of points in a domain by the advancing front technique is an order of magnitude faster than the unstructured mesh for a 3D configuration. The Navier–Stokes solver is based on the artificial compressibility approach and the numerical METHODology is based on the higher-order characteristic-based (CB) discretization. The main objective of this research is to use the CB scheme in order to prevent instabilities. Using this inherent upwind technique for estimating convection variables at the mid-point, no artificial viscosity is required at high Reynolds number. The Taylor LEAST-SQUARES METHOD was used for the calculation of spatial derivatives with normalized Gaussian weight functions. An explicit four stage Runge-Kutta scheme with modified coefficients was used for the discretized equations. To accelerate convergence, local time stepping was used in any explicit iteration for steady state test cases and the residual smoothing techniques were used to converge acceleration. The capabilities of the developed 2D incompressible Navier-Stokes code with the proposed MESHLESS METHOD were demonstrated by flow computations in a lid-driven cavity at four Reynolds numbers. The obtained results using the new proposed scheme indicated a good agreement with the standard benchmark solutions in the literature. It was found that using the third order accuracy for the proposed METHOD could be more efficient than its second order accuracy discretization in terms of computational time.